Selección de sitios para el establecimiento de una red de estaciones meteoclimáticas en República Dominicana, usando AHP y análisis de autocorrelación espacial
José Martínez
Michela Izzo
Introducción
El método de selección de alternativas multicriterio AHP (Analytic Hierarchy Process) se fundamenta en la teoría de la toma de decisiones multicriterio (MCA) y la teoría de la jerarquía analítica. Fue desarrollado por Thomas Saaty en la década de 1970 (Thomas L. Saaty 1977), con varias revisiones posteriores (Thomas L. Saaty 2001; Thomas L. Saaty and Tran 2007), y se utiliza para tomar decisiones cuando se deben considerar múltiples criterios y alternativas. Tradicionalmente, el método AHP se ha utilizado en investigaciones del ámbito de las ingenierías, ciencias sociales, económicas y empresariales, e igualmente en la toma de decisiones donde intervienen datos geoespaciales (Thomas L. Saaty 2013; Darko et al. 2019; Podvezko 2009; Subramanian and Ramanathan 2012; Breaz, Bologa, and Racz 2017). Recientemente, fue usado de forma eficiente en la selección de sitios idóneos para la instalación de estaciones meteoclimáticas en Perú (Rojas Briceño et al. 2021).
El método AHP consiste en descomponer un problema complejo en una estructura jerárquica de criterios y subcriterios, para luego comparar distintas alternativas en función de cada uno de dichos criterios. El proceso se realiza en varias etapas, que incluyen, identificar los objetivos y criterios relevantes para el problema, crear una estructura jerárquica de los criterios y subcriterios, comparar los criterios y subcriterios mediante una matriz de comparación en parejas (paso clave), calcular los valores de prioridad de cada criterio (paso clave), comparar las alternativas, calcular los valores de prioridad de cada alternativa en función de cada criterio y, finalmente, calcular los valores totales de prioridad de cada alternativa.
El método AHP es ampliamente utilizado en la toma de decisiones y en la planificación estratégica, ya que permite elegir entre varias opciones considerando valoraciones de criterios, y porque tiene en cuenta la importancia relativa de los criterios elegidos. Esta importancia relativa se asigna, normalmente, por medio de consultas hechas a personas con experiencia en el área de conocimiento donde se enmarque el problema en cuestión.
En este estudio, aplicamos AHP para seleccionar sitios idóneos donde instalar estaciones meteoclimáticas en República Dominicana, garantizando la eficiencia de la red, maximizando recursos y evitando redundancia información. Para ello, nos apoyamos tanto en fuentes de información geoespacial sistemáticamente producidas, como en consultas a personas con experiencia en temas climáticos y meteorológicos.
Materiales y método
El método AHP se utiliza para seleccionar la mejor opción entre diferentes alternativas, utilizando criterios de selección ponderados por personas con conocimiento del problema (Thomas L. Saaty 2013). Las repuestas originales normalmente deben organizarse y recodificarse y, posteriormente, se debe evaluar su consistencia. A continuación, se seleccionan las respuestas consistentes, o se ajustan las inconsistentes, y se establece la ponderación de criterios. Finalmente, la ponderación definida, se aplica a las fuentes de información disponible para obtener una lista de alternativas, de entre las cuales, se selecciona la más idónea de acuerdo con los criterios definidos.
Tanto el diseño de los formularios, como el procesamiento de respuestas y la ponderación de criterios, los realizamos empleando lenguajes de programación. Para diseñar los formularios, empleamos paquetes y funciones de Python, mientras que para los análisis nos auxiliamos del paquete ahpsurvey y otros del entorno de programación estadística R, diseñado para tales fines (Cho 2019; R Core Team 2021; Wickham et al. 2019). Describimos estos pasos detalladamente en la sección Información suplementaria.
Resultados
Importancia de los criterios
Las tablas 1 y 2 muestran las preferencias individuales y agregadas, respectivamente, de las personas entrevistadas cuyas respuestas fueron consistentes.
kable_prefind <- flujo_completo_ahp$indpref %>%
mutate(`Persona consultada` = cr_indicador[cr_indicador[,1]==1, 'Persona consultada']) %>%
relocate(`Persona consultada`) %>%
estilo_kable(titulo = 'Preferencias individuales',
cubre_anchura = F) %>%
kable_styling(position = 'left') %>%
column_spec(column = 1:2, width = "10em")
kable_prefind| Persona consultada | acce | temp | pluv | habi | agua | pend | inso | elev | CR | top1 | top2 | top3 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0.04 | 0.22 | 0.32 | 0.16 | 0.02 | 0.02 | 0.08 | 0.14 | 0.09 | acce_pend | pluv_elev | pend_inso |
| 2 | 0.09 | 0.23 | 0.25 | 0.07 | 0.07 | 0.03 | 0.11 | 0.16 | 0.07 | habi_agua | habi_inso | pluv_agua |
| 4 | 0.06 | 0.05 | 0.28 | 0.05 | 0.06 | 0.06 | 0.31 | 0.14 | 0.07 | temp_inso | habi_inso | temp_pend |
| 9 | 0.10 | 0.19 | 0.25 | 0.07 | 0.09 | 0.04 | 0.22 | 0.04 | 0.06 | habi_agua | agua_elev | acce_elev |
kable_prefagg <- flujo_completo_ahp$aggpref %>%
as.data.frame() %>%
rownames_to_column('Variable') %>%
mutate(Variable = factor(Variable, labels = variables[sort(names(variables))])) %>%
arrange(desc(AggPref)) %>%
estilo_kable(titulo = 'Preferencias agregadas',
cubre_anchura = F) %>%
kable_styling(position = 'left') %>%
column_spec(column = 1:2, width = "10em")
kable_prefagg| Variable | AggPref | SD.AggPref |
|---|---|---|
| estacionalidad pluviométrica | 0.27 | 0.04 |
| horas de insolación | 0.18 | 0.11 |
| estacionalidad térmica | 0.17 | 0.08 |
| elevación | 0.12 | 0.05 |
| heterogeneidad de hábitat | 0.09 | 0.05 |
| distancia a accesos | 0.07 | 0.03 |
| distancia a cuerpos de agua | 0.06 | 0.03 |
| pendiente | 0.04 | 0.02 |
Mapas de los subcriterios
source('R/funciones.R')
library(sf)
library(kableExtra)
res_h3 <- 7 #Escribir un valor entre 4 y 7, ambos extremos inclusive
ruta_ez_gh <- 'https://raw.githubusercontent.com/geofis/zonal-statistics/'
# ez_ver <- 'da5b4ed7c6b126fce15f8980b7a0b389937f7f35/'
ez_ver <- 'd7f79365168e688f0d78f521e53fbf2da19244ef/'
ind_esp_url <- paste0(ruta_ez_gh, ez_ver, 'out/all_sources_all_variables_res_', res_h3, '.gpkg')
ind_esp_url## [1] "https://raw.githubusercontent.com/geofis/zonal-statistics/d7f79365168e688f0d78f521e53fbf2da19244ef/out/all_sources_all_variables_res_7.gpkg"if(!any(grepl('^ind_esp$', ls()))){
ind_esp <- st_read(ind_esp_url, optional = T, quiet = T)
st_geometry(ind_esp) <- "geometry"
ind_esp <- st_transform(ind_esp, 32619)
}
if(!any(grepl('^pais_url$', ls()))){
pais_url <- paste0(ruta_ez_gh, ez_ver, 'inst/extdata/dr.gpkg')
pais <- invisible(st_read(pais_url, optional = T, layer = 'pais', quiet = T))
st_geometry(pais) <- "geometry"
pais <- st_transform(pais, 32619)
}
if(!any(grepl('^ind_esp_inters$', ls()))){
ind_esp_inters <- st_intersection(pais, ind_esp)
colnames(ind_esp_inters) <- colnames(ind_esp)
ind_esp_inters$area_sq_m <- units::drop_units(st_area(ind_esp_inters))
ind_esp_inters$area_sq_km <- units::drop_units(st_area(ind_esp_inters))/1000000
}
if(!any(grepl('^ind_esp_inters$', ls())) && interactive()){
print(ind_esp_inters)
}Distancia a accesos.
# Objeto que acogerá nombres de objetos
objetos <- character()objeto <- 'osm_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'OSM-DIST mean',
umbrales = c(50, 200, 500, 5000),
nombre = variables[[1]],
ord_cat = 'nin_rev')
)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 12.76 289.04 534.24 1243.20 1412.68 32795.66 2get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| distancia a accesos intervalos | distancia a accesos etiquetas | distancia a accesos puntuación |
|---|---|---|
| [12.8,50] | no idóneo | 1 |
| (50,200] | altamente idóneo | 4 |
| (200,500] | moderadamente idóneo | 3 |
| (500,5e+03] | marginalmente idóneo | 2 |
| (5e+03,3.28e+04] | no idóneo | 1 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Estacionalidad térmica.
objeto <- 'tseasonizzo_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'TSEASON-IZZO mean',
umbrales = c(1.1, 1.3, 1.5),
nombre = variables[[2]],
ord_cat = 'ni')
)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.573 1.195 1.301 1.326 1.477 1.866 105get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| estacionalidad térmica intervalos | estacionalidad térmica etiquetas | estacionalidad térmica puntuación |
|---|---|---|
| [0.573,1.1] | no idóneo | 1 |
| (1.1,1.3] | marginalmente idóneo | 2 |
| (1.3,1.5] | moderadamente idóneo | 3 |
| (1.5,1.87] | altamente idóneo | 4 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Estacionalidad pluviométrica.
objeto <- 'pseasonizzo_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'PSEASON-IZZO mean',
umbrales = c(30, 40, 50),
nombre = variables[[3]],
ord_cat = 'ni')
)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 19.51 31.14 43.51 42.76 52.55 89.60 105get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| estacionalidad pluviométrica intervalos | estacionalidad pluviométrica etiquetas | estacionalidad pluviométrica puntuación |
|---|---|---|
| [19.5,30] | no idóneo | 1 |
| (30,40] | marginalmente idóneo | 2 |
| (40,50] | moderadamente idóneo | 3 |
| (50,89.6] | altamente idóneo | 4 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)
# Para comparar con CHELSA
# objeto <- 'chbio15_rcl'
# assign(
# objeto,
# generar_resumen_grafico_estadistico_criterios(
# variable = 'CH-BIO bio15 precipitation seasonality',
# umbrales = c(300, 400, 500),
# nombre = 'Estacionalidad pluviométrica',
# ord_cat = 'ni')
# )
# get(objeto)[c('violin', 'mapa_con_pais', 'intervalos_y_etiquetas_kable')]
# # clipr::write_clip(get(objeto)$intervalos_y_etiquetas)Heterogeneidad de hábitat.
objeto <- 'hethab_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'GHH coefficient_of_variation_1km',
umbrales = c(300, 450, 600),
nombre = variables[[4]],
ord_cat = 'in')
)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 317.6 399.6 481.4 529.5 3563.4get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| heterogeneidad de hábitat intervalos | heterogeneidad de hábitat etiquetas | heterogeneidad de hábitat puntuación |
|---|---|---|
| [0,300] | altamente idóneo | 4 |
| (300,450] | moderadamente idóneo | 3 |
| (450,600] | marginalmente idóneo | 2 |
| (600,3.56e+03] | no idóneo | 1 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Distancia a cuerpos de agua y humedales.
objeto <- 'wbwdist_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'WBW-DIST mean',
umbrales = c(1000, 2000, 3000),
nombre = variables[[5]],
ord_cat = 'ni')
)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 2698 6069 7134 10545 26424get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| distancia a cuerpos de agua intervalos | distancia a cuerpos de agua etiquetas | distancia a cuerpos de agua puntuación |
|---|---|---|
| [0,1e+03] | no idóneo | 1 |
| (1e+03,2e+03] | marginalmente idóneo | 2 |
| (2e+03,3e+03] | moderadamente idóneo | 3 |
| (3e+03,2.64e+04] | altamente idóneo | 4 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Pendiente.
objeto <- 'slope_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'G90 Slope',
umbrales = c(3, 9, 15),
nombre = variables[[6]],
ord_cat = 'in')
)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.370 4.489 6.773 10.840 32.705get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| pendiente intervalos | pendiente etiquetas | pendiente puntuación |
|---|---|---|
| [0,3] | altamente idóneo | 4 |
| (3,9] | moderadamente idóneo | 3 |
| (9,15] | marginalmente idóneo | 2 |
| (15,32.7] | no idóneo | 1 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Horas de insolación.
objeto <- 'insol_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'YINSOLTIME mean',
umbrales = c(3900, 4100, 4300),
nombre = variables[[7]],
ord_cat = 'ni')
)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 3176 4092 4296 4232 4421 4483 104get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| horas de insolación intervalos | horas de insolación etiquetas | horas de insolación puntuación |
|---|---|---|
| [3.18e+03,3.9e+03] | no idóneo | 1 |
| (3.9e+03,4.1e+03] | marginalmente idóneo | 2 |
| (4.1e+03,4.3e+03] | moderadamente idóneo | 3 |
| (4.3e+03,4.48e+03] | altamente idóneo | 4 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Elevación.
objeto <- 'ele_rcl'
assign(
objeto,
generar_resumen_grafico_estadistico_criterios(
variable = 'CGIAR-ELE mean',
umbrales = c(200, 400, 800),
nombre = variables[[8]],
ord_cat = 'ni')
)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -42.00 57.76 186.78 386.79 542.70 2791.69 35get(objeto)[c('violin', 'mapa_con_pais')]## $violin
##
## $mapa_con_pais
get(objeto)[['intervalos_y_etiquetas_kable']]| elevación intervalos | elevación etiquetas | elevación puntuación |
|---|---|---|
| (200,400] | marginalmente idóneo | 2 |
| (400,800] | moderadamente idóneo | 3 |
| (800,2.79e+03] | altamente idóneo | 4 |
| [-42,200] | no idóneo | 1 |
# clipr::write_clip(get(objeto)$intervalos_y_etiquetas)
if(!objeto %in% objetos) objetos <- c(objetos, objeto)Umbrales de criterios
Los umbrales elegidos para definir las puntuaciones de criterios, están recogidos en la tabla (ver tabla 11).
puntuaciones_umbrales <- map(objetos, function(x) get(x)[['intervalos_y_etiquetas']] %>%
pivot_longer(cols = -matches('puntuación|etiquetas'), names_to = 'criterio') %>%
mutate(criterio = gsub(' intervalos', '', criterio)) %>%
group_by(across(all_of(matches('etiquetas|criterio')))) %>%
summarise(value = paste(value, collapse = ' y ')) %>%
pivot_wider(names_from = contains('etiquetas'), values_from = value) %>%
select(criterio, `altamente idóneo`, `moderadamente idóneo`, `marginalmente idóneo`, `no idóneo`)
) %>% bind_rows()
readODS::write_ods(puntuaciones_umbrales, 'fuentes/umbrales-criterios-ahp/puntuaciones.ods')
puntuaciones_umbrales %>% kable(format = 'html', escape = F, booktabs = T, digits = 2,
caption = 'Puntuaciones de criterios para la selección de sitios de estaciones meteoclimáticas') %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = T)| criterio | altamente idóneo | moderadamente idóneo | marginalmente idóneo | no idóneo |
|---|---|---|---|---|
| distancia a accesos | (50,200] | (200,500] | (500,5e+03] | [12.8,50] y (5e+03,3.28e+04] |
| estacionalidad térmica | (1.5,1.87] | (1.3,1.5] | (1.1,1.3] | [0.573,1.1] |
| estacionalidad pluviométrica | (50,89.6] | (40,50] | (30,40] | [19.5,30] |
| heterogeneidad de hábitat | [0,300] | (300,450] | (450,600] | (600,3.56e+03] |
| distancia a cuerpos de agua | (3e+03,2.64e+04] | (2e+03,3e+03] | (1e+03,2e+03] | [0,1e+03] |
| pendiente | [0,3] | (3,9] | (9,15] | (15,32.7] |
| horas de insolación | (4.3e+03,4.48e+03] | (4.1e+03,4.3e+03] | (3.9e+03,4.1e+03] | [3.18e+03,3.9e+03] |
| elevación | (800,2.79e+03] | (400,800] | (200,400] | [-42,200] |
Representación reclasificada de criterios y puntuaciones agregadas
Unir los vectoriales de cada criterio y representar mapa.
all_criteria <- map(objetos[2:length(objetos)], ~ get(.x)[['vectorial']] %>% st_drop_geometry) %>%
prepend(list(get(objetos[1])[['vectorial']])) %>%
reduce(left_join, by = "hex_id")
all_criteria %>% st_write('out/intervalos_etiquetas_puntuaciones_AHP_criterios_separados.gpkg', delete_dsn = T)## Deleting source `out/intervalos_etiquetas_puntuaciones_AHP_criterios_separados.gpkg' using driver `GPKG'
## Writing layer `intervalos_etiquetas_puntuaciones_AHP_criterios_separados' to data source
## `out/intervalos_etiquetas_puntuaciones_AHP_criterios_separados.gpkg' using driver `GPKG'
## Writing 13152 features with 25 fields and geometry type Unknown (any).Mapas puntuaciones reclasificadas de cada criterio.
paleta <- c("altamente idóneo" = "#018571", "moderadamente idóneo" = "#80cdc1",
"marginalmente idóneo" = "#dfd2b3", "no idóneo" = "#a6611a")
all_criteria_mapa <- all_criteria %>%
select(all_of(contains('etiquetas'))) %>%
rename_with(~ stringr::str_replace(.x,
pattern = ' etiquetas',
replacement = ''),
matches('etiquetas')) %>%
pivot_longer(cols = -geometry) %>%
ggplot +
aes(fill = value) +
geom_sf(lwd=0) +
scale_fill_manual(values = paleta) +
labs(title = paste('Reclasificación de valores de criterios')) +
geom_sf(data = pais, fill = 'transparent', lwd = 0.5, color = 'grey50') +
facet_wrap(~ name, ncol = 2) +
theme_bw() +
theme(
legend.position = 'bottom',
legend.key.size = unit(0.5, 'cm'), #change legend key size
legend.key.height = unit(0.5, 'cm'), #change legend key height
legend.key.width = unit(0.5, 'cm'), #change legend key width
legend.title = element_blank(), #change legend title font size
legend.text = element_text(size=2) #change legend text font size
)
if(interactive()) dev.new()
all_criteria_mapa